function OO_minloss = cf_minloss(low_multiplier, high_multiplier, ...
    cost2_0, cost2_1_all, D_input, D, N,...
    D_enroll_st, ...
    alpha, ...
    Delta_xi_all, cost1_e_all, mu_f_all, sigma_f_all,...
    K, pr_y, N_y, resource_y1, resource_y2, y_medi, x_fc, pr_fc, ...
    rho, crra, sigma_c1, beta_c, beta_f, T1, T2, B1_c, B2_c, B1_f, B2_f, delta, ...
    lr_target, N_f, entry_exitflag_vec, id_ist_old, myopic, Er, cost_weight, use_cost_weight, root_dir)

% loop over different regimes
for ss=1:3
    tic0=tic;
    disp(' ')

    if ss==1 % low mins loss reg costs
        disp('***************************')
        disp('ss=1, Initial rate regulatory cost small')

        cost1_e_cf = low_multiplier*cost1_e_all;

    elseif ss==2 % high min loss reg costs
        disp('***************************')
        disp('ss=2, Initial rate regulatory cost large')

        cost1_e_cf = high_multiplier*cost1_e_all;

    elseif ss==3 % baseline
        disp('***************************')
        disp('ss=3, Baseline regime')

        cost1_e_cf = cost1_e_all;

    end

    % initial guess on eqm outcomes: use observed data
    p1_ini = D_input.p1; % NN x 1
    p2_ini = D_input.p2; % NN x KM
    n_ini_unmerge = D_input.n_insurers; % NN x 1: 1 for majors, # fringes for the fringe
    n_ini = zeros(N.st,1); % N_st x 1. bring down to market level
    for mm=1:N.st
        n_ini(mm) = unique(n_ini_unmerge(D_input.mkt==mm & D_input.major==0));
    end

    % solve the counterfactual eqm
    [p1_cf, p2_cf, n_cf, ...
        profit_cf, non_converge_all_cf, solve_market_all_cf, fsolve_all_cf, ...
        EV_y_cf, EV_cf, enrollment_cf, enrollment_m_cf, enrollment_f_cf, enrollment_per_m_cf, enrollment_per_f_cf, non_converge_cf, solve_market_cf,...
        s1_mkt_cf, s2_mkt_cf, rs_share_cf]...
        = solve_eq(alpha, ...
        Delta_xi_all, cost1_e_cf, cost2_1_all,  mu_f_all, sigma_f_all,...
        cost2_0, D_input, N, ...
        K, pr_y, N_y, resource_y1, resource_y2, y_medi, x_fc, pr_fc, ...
        rho, crra, sigma_c1, beta_c, beta_f, T1, B1_c, B2_c, B1_f, B2_f, delta, ...
        lr_target, N_f, entry_exitflag_vec, id_ist_old, myopic, Er,...
        p1_ini, p2_ini, n_ini, cost_weight, use_cost_weight);

    % save all outcomes
    OO_minloss(ss).p1_cf=p1_cf;
    OO_minloss(ss).p2_cf=p2_cf;
    OO_minloss(ss).n_cf=n_cf;
    OO_minloss(ss).profit_cf=profit_cf;
    OO_minloss(ss).non_converge_all_cf = non_converge_all_cf;
    OO_minloss(ss).solve_market_all_cf = solve_market_all_cf;
    OO_minloss(ss).fsolve_all_cf = fsolve_all_cf;
    OO_minloss(ss).EV_y_cf = EV_y_cf;
    OO_minloss(ss).EV_cf = EV_cf;
    OO_minloss(ss).enrollment_cf = enrollment_cf;
    OO_minloss(ss).enrollment_m_cf = enrollment_m_cf;
    OO_minloss(ss).enrollment_f_cf = enrollment_f_cf;
    OO_minloss(ss).enrollment_per_m_cf = enrollment_per_m_cf;
    OO_minloss(ss).enrollment_per_f_cf = enrollment_per_f_cf;
    OO_minloss(ss).non_converge_cf=non_converge_cf;
    OO_minloss(ss).solve_market_cf=solve_market_cf;
    OO_minloss(ss).s1_mkt_cf = s1_mkt_cf;
    OO_minloss(ss).s2_mkt_cf = s2_mkt_cf;
    OO_minloss(ss).rs_share_cf = rs_share_cf;

    toc0 = toc(tic0);
    disp(' ')
    disp(['Time elapsed for current ss (min) = ', num2str(toc0/(60), 5)]);
end

% save
save('cf_minloss.mat')

